Table of Contents
- Why build phylogenies?
- How phylogenies are reconstructed
- Some terminology
- Simulating evolution
- Visualizing trees with ete3
- Distance-based approaches to phylogenetic reconstruction
- Distances and distance matrices
- Alignment-free distances between sequences
- Alignment-based distances between sequences
- Jukes-Cantor correction of observed distances between sequences
- Phylogenetic reconstruction with UPGMA
- Phylogenetic reconstruction with neighbor-joining
- Limitations of distance-based approaches
- Bootstrap analysis
- Parsimony-based approaches to phylogenetic reconstruction
- Statistical approaches to phylogenetic reconstruction
- Rooted versus unrooted trees
In this chapter we'll begin to explore the goals, approaches, and challenges for creating phylogenetic trees, or phylogenies. Phylogenies, such as the two presented in Figure 1, represent hypotheses about the evolutionary history of a groups, who are represented by the tips in the tree. You can explore an interactive version of the three-domain tree presented in Figure 1b online, through the Interactive Tree of Life project.
- https://en.wikipedia.org/wiki/Computational_phylogenetics
- https://en.wikipedia.org/wiki/Distance_matrices_in_phylogeny
- https://en.wikipedia.org/wiki/Phylogenetic_tree
- https://en.wikipedia.org/wiki/Tree_of_life_(biology)
- https://en.wikipedia.org/wiki/Cladistics
- https://en.wikipedia.org/wiki/Systematics
- https://en.wikipedia.org/wiki/Dendrogram
- https://en.wikipedia.org/wiki/Cladogram
- https://en.wikipedia.org/wiki/Phylogenetic_nomenclature
- https://en.wikipedia.org/wiki/Phylogenetic_network
|
|
| Figure 1a: Evolutionary tree presented by Charles Darwin in On the Origin of Species. Figure 1b: A hypothesis of evolutionary relationships between the three domains of life. This image was created by the Norman Pace Laboratory. | |
Reconstructing the phylogeny of a group of individuals is useful for many reasons. Probably the most obvious of these is understanding the evolutionary relationship between a group of organisms. For example, over the past half-century we've gained great insight into the evolution of our species, Homo sapiens, by studying features of our closest relatives, both extant (still existing) organisms, such as the Pan (chimpanzees) and Gorilla genera, and extinct (no longer living) species, including Homo neanderthalensis, Homo erectus, Homo habilus, and many species in the Australopithecus genus. (The Smithsonian Museum's Human Origins Initiative is an excellent resource for learning more about the fascinating subject of human evolution.)
We've also improved our understanding of the deeper branches in the tree of life.¶
For example, Woese and Fox (1977) used phylogenetic reconstruction to first illustrate that the "prokaryotes" really represented two ancient lineages which they called the eubacteria and the archaebacteria, which ultimately led to the proposal of a "three domain" tree of life composed of the three deep branching linages, the archaea, the bacteria, and the eucarya (Woese, Kandler and Wheelis (1990)).
Phylogenetic trees such as these are also useful for understanding evolution itself.
In fact, they're so useful that the single image that Charles Darwin found important enough to include in On the Origin of Species was the phylogenetic tree presented in Figure 1a.¶
The individuals represented at the tips of our trees don't necessarily have to be organisms though. In another important application of phylogenetic trees, we can study the evolution of genes, which can help us gain a deeper understanding of gene function. Specifically, we can learn about families of related genes.
A classic example of this is the globin family, which includes the proteins hemoglobin and myoglobin, molecules that can reversibly bind oxygen (meaning they can bind to it, and then let go of it). You've probably heard of hemoglobin (if not globins in general), as this molecule binds to oxygen where it is present in high concentration (such as in your lung) and releases it where it is present in low concentration (such as in the bicep, where it is ultimately used to power your arm). Hemoglobin and myoglobin are paralogs, meaning that they are related by a gene duplication and subsequent divergence. If you were to compare an unknown globin sequence to either of these you could detect homology, but a tree such as the one present in Figure 2, would help you understand the type of homologous relationship (i.e., whether it was orthology or paralogy).
Phylogenetic trees are used for many other diverse applications in bioinformatics, so it's therefore important that a bioinformatician have an understanding of they are built and how they should be interpreted.¶
An additional application that we'll cover in this text is comparing the composition of communities of organisms, but we'll come back to that later.
Phylogenies are reconstructed using a variety of different algorithms, some of which we'll cover in this chapter. These algorithms all work by comparing a set of features of organisms, and inferring the evolutionary distance between those organisms based on the similarity of their features. The features that are compared can be nearly anything that is observable, either from extant organisms or fossilized representatives of extinct organisms.
As an example, let's consider the reconstruction of the phylogeny of spiders (the order Araneae), a hypothesis of which is presented in Figure 3. Of the extant spiders, some are orb-weavers (meaning they spin circular, sticky webs), and others are not. Entomologists have debated whether orb-weaving is a monophyletic trait (meaning that it evolved one time), or whether it is polyphyletic (meaning that it evolved multiple times, such as flight, which has evolved independently in birds, flying dinosaurs, insects, and mammals). If orb-weaving is monophyletic, it would mean that over the course of evolution, extant spiders which don't weave orb webs have lost that ability. Some researchers doubt this as it's a very effective means of catching prey, and losing that ability would likely constitute an evolutionary disadvantage. If orb-weaving is polyphyletic, it would means that in at least two different spider lineages, this trait arose independently, which other researchers consider to be very unlikely due to the complexity of engineering these webs. Examples of the evolution of monophyletic and polyphyletic traits are presented in Figures 3 and 4, respectively.
Kahoot ------------¶
Earlier work on understanding the relations between the spider lineages focused on comparing traits that entomologists would observe, for example by watching spiders in action or by dissecting them. For example, in 1986 through 1991, Johnathan Coddington published several studies that tabulated and compared about sixty features of 32 spider taxa (Coddington J. 1986. The monophyletic origin of the orb web. In: Shear W, ed. Spiders: webs, behavior, and evolution. Stanford, California: Stanford University Press. 319-363; Coddington JA. 1991. Cladistics and spider classification: araneomorph phylogeny and the monophyly of orbweavers (Araneae: Araneomorphae; Orbiculariae) Acta Zoologica Fennica 190:75-87). Features included whether a spider wraps its prey when it attacks (a behavioral trait), and how branched the spider's trachea is (a morphological trait). By determining which spiders were more similar and different across these traits, Dr. Coddington provided early evidence for the hypothesis that orb-weaving is an ancient monophyletic trait.
More recently, several research teams have used features of spider genomes to reconstruct the spider phylogeny (Bond et al., 2014, Garrison et al., 2016).¶
- Using this approach, the features become the nucleotides observed at particular positions in the genome, which are observed first by sequencing specific genes that the researchers target that are present in all members of the group, and then aligning those sequences with multiple sequence alignment.
This has several advantages over feature matrices derived from morphological and behavioral traits, including that many more features can be observed.¶
For example, (Garrison et al., 2016), compared approximately 700,000 amino acid positions from nearly 4000 loci around the genomes of 70 spider taxa. Compare the number of features here to the number mentioned in the previous paragraph. These phylogenomic studies have further supported the idea that orb-weaving is an ancient monophyletic trait, and have provided much finer scale information on the evolution of spiders. Supported by these data, researchers hypothesize that the loss of orb-weaving might not be that surprising. While it does provide an effective means of catching flying insects, many insects which are potential prey for spiders don't fly. Further, orb webs may attract predators of spiders, as they are easily observable signals of where a spider can be found.
For the remainder of this chapter, we'll consider methods for phylogenetic reconstruction that use genome sequence data as features.
Terminal nodes or tips typically represent extant organisms, also frequently called operational taxonomic units or OTUs.¶
- OTU is a generic way of referring to a grouping of organisms (such as a species, a genus, or a phylum), without specifically identifying what that grouping is.
Internal nodes in a phylogenetic tree represent hypothetical ancestors.¶
- We postulate their existence but often don't have direct evidence.
- The root node is the internal node from which all other nodes in the tree descend. This is often referred to as the last common ancestor (LCA) of the OTUs represented in the tree. In a universal tree of life, the LCA is often referred to as LUCA, the last universal common ancestor. All nodes in the tree can be referred to as OTUs.
Branches connect the nodes in the tree, and generally represent time or some amount of evolutionary change between the OTUs.¶
- The specific meaning of the branches will depend on the method that was used to build the phylogenetic tree.
A clade in a tree refers to some node (either internal or terminal) and all nodes descending from it (i.e., moving away from the root toward the tips).¶
We'll use all of these terms below as we begin to explore phylogenetic trees.
Table of Contents
Before we jump into how to reconstruct a phylogeny from DNA sequence data, we're going to perform a simulation of the process of evolution of a DNA sequence. In this simulation, we're going to model sequence evolution with a Python function, and then we're going to run that function to simulate multiple generations of evolution.
Bioinformatics developers often use simulations to understand how their algorithms work, as they uniquely provide an opportunity to know what the correct answer is. This provides a way to compare algorithms to each other to figure out which performs best under which circumstances. In our simulation we're going to have control over the starting sequence, and the probability of incurring a substitution mutation or an insertion/deletion mutation at each position of the sequence in each generation.
This would, for example, let us understand whether different algorithms for phylogenetic reconstruction are better or worse for more closely or distantly related sequences.¶
Our simulation will works as follows.
- We'll have one function that we primarily interact with called
evolve_generations.- This will take a starting sequence and the number of generations that we want to simulate.
- It will also take the probability that we want a substitution mutation to occur at each position in the starting sequence, and the probability that we want either an insertion or deletion mutation to occur at each position.
- In each generation, every sequence will spawn two new sequences, randomly incurring mutations at the prescribed rates.
- This effectively simulates a clonal process of reproduction, a form of asexual reproduction common in single cellular organisms, where a parent cell divides into two cells, each containing a copy of the parent's genome with some errors introduced.
Let's inspect this code and then run our simulation beginning with a random sequence.
%pylab inline
import numpy as np
import seaborn as sns
import random
#from IPython.core import page
#page.page = print
First we'll look at the function used to simulate the evolution of a single sequence. This is where most of the important evolutionary modeling of sequence evolution happens.
from iab.algorithms import evolve_sequence
%psource evolve_sequence
Next, take a look at the function that models a single generation of a single sequence. This is where the clonal reproduction (i.e., one parent sequence becoming two child sequences) occurs.
from iab.algorithms import evolve_generation
%psource evolve_generation
Finally, take a look at our entry point function. This is where we provide the parameters of our simulation, including the starting sequence, the number of generations, and the mutation probabilities. Notice how each of these functions builds on the prior functions.
from iab.algorithms import evolve_generations
%psource evolve_generations
Now we'll run our simulation. We'll start with a random DNA sequence, and then evolve three generations.
Before running this, can you predict how many child sequences we'll end up with after three generations?¶
When we call evolve_generations, we'll pass the parameter verbose=True. This will tell the function to print out some information throughout the process. This will let us inspect our evolutionary process: something that is impossible to do with real sequences from history, though we can observe microbial evolution in the lab.
from iab.algorithms import random_sequence
import skbio
sequence = random_sequence(skbio.DNA, 50)
sequences = evolve_generations(sequence, generations=3,
substitution_probability=0.1,
indel_probability=0.05,
increased_rate_probability=0.1,
verbose=True)
We now have a new variable, sequences, which contains the child sequences from the last generation.
Take a minute to look at the ids of the parent and child sequences above, and the ids of a couple of the final generation sequences.¶
- These ids are constructed so that each sequence contains the identifiers of its ancestral sequences, and then either
1or2. - Notice that all sequence identifiers start with
0, the identifier of the last common ancestor (or our starting sequence) of all of the sequences.
These identifiers will help us interpret whether the phylogenies that we reconstruct accurately represent the evolutionary relationships between the sequences.
Also, notice that at this point we only have the sequences from the last generation. We no longer have the ancestral sequences (which would correspond to the internal nodes in the tree). This models the real world, were we only have sequences from extant organisms, but not their ancestors.¶
Take a minute to compare the two sequences below. What types of mutations happened over the course of their evolution?
print(len(sequences))
sequences[0]
sequences[-1]
In our simulation, each sequence is directly derived from exactly one sequence from the previous generation.¶
- The evolution of all of the sequences traces back to starting sequence that we provided.
- This means that our final sequences are all homologous.
And because we have modeled this process, we know where each sequence fits in relation to all of the other sequences in the phylogeny.¶
- Our goal with the algorithms we'll study for the rest of this chapter is to reconstruct that phylogeny given only the last generation of sequences.
- We'll use the fact that we know the true phylogeny to help us evaluate the relative performance of the different methods.
Let's simulate 10 generations of sequences here, and then randomly select some of those sequences to work with in the remaining sections of this chapter. For the sake of runtime, I'm going to set our indel_probability=0.0. This means that we won't need to align our selected sequences before constructing trees from them (because with only substitution mutations occurring, homologous positions always remain aligned). Multiple sequence alignment is currently a runtime bottleneck in IAB, so this will let us run this notebook much faster. If you'd like to model insertion/deletions, you can increase the indel_probability, say to 0.005. If you do that, the sequences will be aligned for you in the next cell, but it may take around 30 minutes to run.
from skbio.alignment import global_pairwise_align_nucleotide
from iab.algorithms import progressive_msa
from functools import partial
indel_probability = 0.0
sequences = evolve_generations(sequence,
generations=10,
substitution_probability=0.03,
indel_probability=0.0,
increased_rate_probability=0.1,
verbose=False)
sequences = random.sample(sequences, 25)
if indel_probability == 0:
sequences_aligned = sequences
else:
gpa = partial(global_pairwise_align_nucleotide,
penalize_terminal_gaps=True)
sequences_aligned = progressive_msa(sequences,
pairwise_aligner=gpa)
While simulations are extremely powerful for comparing algorithms, they can also be misleading. This is in part because when we model evolution we simplify the evolutionary process, and the model is stochastic. One assumption that our simulation is making is that "bursts" of evolution (i.e., where our substitution rate temporarily increases) are restricted to only a single generation. A child is no more or less like to have an increased rate of substitutions if its parent did. This may or may not be an accurate model. Imagine in the real-world that the environment changed drastically for some descendants (for example, if a geological event created new thermal vents in a lake they inhabited, resulting in an increase in mean water temperature), but not for others. The lineages who experience the environmental change might have an increased rate of substitutions as their genomes adapt to the new environment. The increased substitution rate may persist for multiple generations, or it may not. This is one limitation of our simulation that we know about, but because we don't have a perfect understanding of sequence evolution, there are limitations that we don't know about.
When using simulations, it's important to understand what assumptions the simulation makes so you know what it can tell you about and what it can't tell you about. You'll want to consider that when determining how confident you are in the results of an evaluation based on simulation.
Why question: What are some other assumptions that are being made by the evolutionary simulation presented here? There are many, so take a minute to list a few.¶
On the opposite end of the spectrum from simulations for algorithm comparison is comparisons based on real data, where we don't know groud-truth.
- The trade-off is that, while there are no assumptions being made in real data with respect to the composition of the sequences, with real data we don't usallly know what the right answer (in our case, the correct phylogeny) is (unless we record it in the a microbiology lab), so it's harder to determine which algorithms are doing better or worse.
The take-away message here is that neither approach is perfect, and often researchers will use a combination of simulated and real data to evaluate algorithms.
As we now start computing phylogenetic trees, we're going to need a way to visualize them. We'll use the ete3 Python package for this, and in the next cell we'll configure the TreeStyle which is used to define how the trees we visualize will look. If you'd like to experiment with other views, you can modify the code in this cell according to the ete3 documentation. If you come up with a nicer style, I'd be interested in seeing that in a pull request. You can post screenshots to IAB issue #213 before creating a pull request so I can see what the new style looks like.
import ete3
help(ete3.Tree)
ts = ete3.TreeStyle()
ts.show_leaf_name = True
ts.scale = 250
ts.branch_vertical_margin = 15
We can apply this TreeStyle to a random tree as follows. Any changes that you make to the TreeStyle above will impact this tree and all following trees in this chapter. Experiment with this - it's fun!
# Make a tree object
t = ete3.Tree()
# Generates a random tree of 10
t.populate(10)
# Plot the tree
t.render("%%inline", tree_style=ts)
Table of Contents
- Distances and distance matrices
- Alignment-free distances between sequences
- Alignment-based distances between sequences
- Jukes-Cantor correction of observed distances between sequences
- Phylogenetic reconstruction with UPGMA
- Phylogenetic reconstruction with neighbor-joining
- Limitations of distance-based approaches
The next approaches we'll take for phylogenetic reconstruction rely on computing distances between sequences. We've previously discussed distances between sequences in a few places in the text. We'll begin this section by formalizing the term distance, and introducing the concept of a distance matrix.
Generally, we can consider a distance between a pair of objects is a measure of their dissimilarity.¶
- There isn't one single definition of the distance between two objects.
- For example, if your two objects are the cities Flagstaff, Arizona and Boulder, Colorado you might measure the distance between them as length of shortest line that connects them. This would be the Euclidean distance between the two cities.
- If you're trying to travel from one city to another however, that might not be the most relevant distance measure. Instead, you might be interested in something like their Manhatten distance, which in this would be more similar to a measure of the shortest route that you could travel between the two cities using the Interstate Highway system. You can see how this is different than the shortest line connecting the two cities.
An illustration comparing the taxicab metric to the Euclidean metric on the plane:
- According to the taxicab metric all three pictured paths (red, yellow, and blue) have the same length (12).
- According to the Euclidean metric, the green path has length and the unique shortest path is $$6\sqrt{2} = 8.49$$
Similarly, different distance metrics will be relevant and not relevant for different types of objects. Clearly the way you measure distances between cities is very different than the way we measure distances between biological sequences. However, the underlying concept of a distance between two objects is the same.¶
https://en.wikipedia.org/wiki/Metric_(mathematics)
https://en.wikipedia.org/wiki/Metric_space
- In mathematics, a metric or distance function is a function that defines a distance between each pair of elements of a set.
- A set with a metric is called a metric space.
https://en.wikipedia.org/wiki/Similarity_measure
- In statistics and related fields, a similarity measure or similarity function is a real-valued function that quantifies the similarity between two objects.
- Although no single definition of a similarity measure exists, usually such measures are in some sense the inverse of distance metrics:
- they take on large values for similar objects and either zero or a negative value for very dissimilar objects.
Formally, a measure of dissimilarity $d$ between two objects $x$ and $y$ is a distance if it meets these four criteria for all $x$ and $y$:
- $d(x,y) \geq 0$ (non-negativity)
- $d(x,y) = 0\ if\ and\ only\ if\ x = y$ (identity of indiscernibles)
- $d(x,y) = d(y,x)$ (symmetry)
- $d(x,z) \leq d(x,y) + d(y,z)$ (triangle inequality)
When we compute the distances between some number of objects $n$, we'll commonly represent those values in a distance matrix which contains all of those values.
Note the tree on the axes¶
These distance matrices are so common in bioinformatics that scikit-bio defines a DistanceMatrix object that provides a convenient interface for working with these data.
from skbio import DistanceMatrix
from skbio.stats.distance._base import DissimilarityMatrix
#help(DistanceMatrix)
help(DissimilarityMatrix)
We can create one of these objects as follows:
dm = DistanceMatrix([[0.0, 1.0, 2.0],
[1.0, 0.0, 3.0],
[2.0, 3.0, 0.0]],
ids=['a', 'b', 'c'])
We can then access the values in the distance matrix directly, view the distance matrix as a heatmap, and do many other things that facilitate analyzing distances between objects.
print(dm['a', 'b'])
print(dm['b', 'c'])
help(dm.plot)
_ = dm.plot(cmap='Greens')
The conditions of a distance metric listed above lead to a few specific features of distance matrices:¶
- they're symmetric (if you flip the upper triangle over the diagonal, the values are the same as those in the lower triangle),
- hollow (the diagonal is all zeros), and
- all values are greater than or equal to zero.
Which of the conditions listed above results in each of these features of distance matrices?
We've now looked at several ways of computing distances between sequences, some of which have required that the positions in the sequences are directly comparable to one another (i.e., that our sequences are aligned), and some of which haven't.
One alignment-free distance between sequences is the k-mer distance that we worked with in Sequence Homology Searching.¶
We can use the kmer_distance function with scikit-bio as follows to create an skbio.DistanceMatrix object.
from iab.algorithms import kmer_distance
%psource kmer_distance
help(DistanceMatrix.from_iterable)
print(sequences)
kmer_dm = DistanceMatrix.from_iterable(sequences,
metric=kmer_distance,
key='id')
print(kmer_dm)
_ = kmer_dm.plot(cmap='Greens',
title='3mer distances between sequences')
One alignment-based distance metric that we've looked at is Hamming distance. This would be considered an alignment-based approach because is does consider the order of the characters in the sequence by comparing a character at a position in one sequence only to the character at the corresponding position in the other sequence. We could compute these distances as follows, after first aligning our sequences.
Note: here, we specified a 0 probability of indel mutations, with sequences of identical length, and thus they are already aligned:
from skbio.sequence.distance import hamming
hamming_dm = DistanceMatrix.from_iterable(sequences_aligned,
metric=hamming,
key='id')
_ = hamming_dm.plot(cmap='Greens',
title='Hamming distances between sequences')
Why question: How long should branches be? Does it matter for our conclusions?¶
UPGMA is "ultrametric", meaning that all the terminal nodes (i.e. the sequences/taxa) are equally distance from the root. In molecular terms, this means that UPGMA assumes a molecular clock, i.e. all lineages are evolving at a constant rate. In practical terms, this means that you can construct a distance scale bar and all the terminal nodes will be level at position 0.0, representing the present. In this example, the scale bar is shown on the right-hand side.
The Hamming distance between aligned sequences, as described above, is simple to calculate, but it is often an underestimate of the actual amount of mutation that has occurred in a sequence. Here's why:¶
- imagine that in one generation $g$, position $p$ of sequence $S1$ undergoes a substitution mutation from
AtoC. - Then, in the next generation $g + 1$, the same position $p$ of sequence $S1$ undergoes a substitution from
CtoT. - Because we can only inspect the modern-day sequences, not their ancestors, it looks like position $p$ has had a single substitution event.
- Similarly, if in generation $g + 1$ position $p$ changed from
Cback toA(a back substitution), we would observe zero substitution mutations at that position even though two had occurred.
To correct for this, the Jukes-Cantor correction is typically applied to the Hamming distances between the sequences.¶
- Where $p$ is the Hamming distance, the corrected genetic distance is computed as $d = -\frac{3}{4} \ln(1 - \frac{4}{3}p)$.
- The derivation of this formula is beyond the scope of this text (you can find it in Inferring Phylogeny by Felsenstein), but it is based on the Jukes-Cantor (JC69) nucleotide substitution model.
The Python implementation of this correction looks like the following. We can apply this to a number of input distance values to understand how it transforms our Hamming distances.
def jc_correction(p):
return (-3 / 4) * np.log(1 - (4 * p / 3))
distances = np.arange(0, 0.70, 0.05)
jc_corrected_distances = list(map(jc_correction, distances))
ax = sns.pointplot(distances, jc_corrected_distances)
ax.set_xlabel('Hamming distance')
ax.set_ylabel('JC-corrected distance')
ax.set_xlim(0)
ax.set_ylim(0)
ax
We can then apply this to a full distance matrix as follows (we'll then print the first row of each).
def jc_correct_dm(dm):
result = np.zeros(dm.shape)
for i in range(dm.shape[0]):
for j in range(i):
result[i, j] = result[j, i] = jc_correction(dm[i, j])
return skbio.DistanceMatrix(result, ids=dm.ids)
jc_corrected_hamming_dm = jc_correct_dm(hamming_dm)
print(hamming_dm[0])
print(jc_corrected_hamming_dm[0])
# Pre-correction:
_ = hamming_dm.plot(cmap='Greens', title='Hamming distances between sequences')
# Post-correction:
_ = jc_corrected_hamming_dm.plot(cmap='Greens', title='JC-corrected Hamming distances between sequences')
Table of Contents
The first algorithm we'll look at for reconstructing phylogenetic trees is called UPGMA, which stands for Unweighted Pair-Group Method with Arithmetic mean. While that name sounds complex, it's actually a straightforward algorithm, which is why we're starting with it. After we work through the algorithm, we'll come back to the name as it'll make more sense then.
UPGMA is a generic hierarchical clustering algorithm. It's not specific to reconstructing biological trees, but rather is used for interpreting any type of distance matrix. It is fairly widely used for building phylogenetic trees, though it's application in phylogenetics is usually restricted to building preliminary trees to "guide" the process of multiple sequence alignment. The reason for this is that it's fast, but it makes some assumptions that don't work well for inferring relationships between organisms, which we'll discuss after working through the algorithm.
UPGMA starts with a distance matrix, and works through the following steps to create a tree.
Step 1: Find the smallest non-zero distance in the matrix and define a clade containing only those members. Draw that clade, and set the total length of the branch connecting the tips to the distance between the tips. The distance between each tip and the node connecting them should be half of the distance between the tips.
Step 2: Create a new distance matrix with an entry representing the new clade created in step 1.
Step 3: Calculate the distance matrix entries for the new clade as the mean distance from each of the tips of the new clade to all other tips in the original distance matrix.
Step 4: If there is only one distance (below or above the diagonal) in the distance matrix, use it to connect the remaining unconnected clades, and stop. Otherwise repeat step 1.
Let's work through these steps for a simple distance matrix representing the distances between five sequences.
_data = np.array([[ 0., 4., 2., 5., 6.],
[ 4., 0., 3., 6., 5.],
[ 2., 3., 0., 3., 4.],
[ 5., 6., 3., 0., 1.],
[ 6., 5., 4., 1., 0.]])
_ids = ['s1', 's2', 's3', 's4', 's5']
master_upgma_dm = skbio.DistanceMatrix(_data, _ids)
print(master_upgma_dm)
Iteration 1¶
Step 1: The smallest non-zero distance in the above matrix is between s4 and s5. So, we'll draw that clade and set each branch length to half of the distance between them.
Step 2: Next, we'll create a new, smaller distance matrix where the sequences s4 and s5 are now represented by a single clade which we'll call (s4, s5). This notation indicates that the corresponding distances are to both s4 and s5.
iter1_ids = ['s1', 's2', 's3', '(s4, s5)']
iter1_dm = [[0.0, 4.0, 2.0, None],
[4.0, 0.0, 3.0, None],
[2.0, 3.0, 0.0, None],
[None, None, None, None]]
Step 3: We'll now fill in the values from the new clade to each of the existing sequences (or clades). The distance will be the mean between each pre-existing clade, and each of the sequences in the new clade. For example, the distance between s1 and (s4, s5) is the mean of the distance between s1 and s4 and s1 and s5:
import numpy as np
# Note: we are not using the iter1_dm, but the original matrix (master_upgma_dm)
s1_s4s5 = np.mean([master_upgma_dm['s1', 's4'], master_upgma_dm['s1', 's5']])
print(s1_s4s5)
Similarly, the distance between s2 and (s4, s5) is the mean of the distance between s2 and s4 and s2 and s5:
s2_s4s5 = np.mean([master_upgma_dm['s2', 's4'], master_upgma_dm['s2', 's5']])
print(s2_s4s5)
And finally, the distance between s3 and (s4, s5) is the mean of the distance between s3 and s4 and the distance between s3 and s5:
s3_s4s5 = np.mean([master_upgma_dm['s3', 's4'], master_upgma_dm['s3', 's5']])
print(s3_s4s5)
We can fill these values in to our iteration 1 distance matrix. Why do we only need to compute three values to fill in seven cells in this distance matrix?
iter1_dm = [[0.0, 4.0, 2.0, s1_s4s5],
[4.0, 0.0, 3.0, s2_s4s5],
[2.0, 3.0, 0.0, s3_s4s5],
[s1_s4s5, s2_s4s5, s3_s4s5, 0.0]]
iter1_dm = DistanceMatrix(iter1_dm, iter1_ids)
print(iter1_dm)
Step 4: Because there is still more than one value below the diagonal in our new distance matrix, we start a new iteration by going back to Step 1 and repeating this process.
Iteration 2¶
Step 1: The smallest non-zero distance in the iteration 1 distance matrix is between s1 and s3. So, we'll draw that clade and set each branch length to half of that distance.
Step 2: We next create a new, smaller distance matrix where the sequences s1 and s3 are now represented by a single clade, (s1, s3).
iter2_ids = ['(s1, s3)', 's2', '(s4, s5)']
iter2_dm = [[None, None, None],
[None, 0.0, 5.5],
[None, 5.5, 0.0]]
Step 3: We'll now fill in the values from the new clade to each of the existing sequences (or clades). Notice that the distance between our new clade and s2 is the mean of two values, but the distance between our new clade and the clade defined in iteration 1 is the mean of four values. Why is that?
# For combining on s2 (index 1)
s2_s1s3 = np.mean([master_upgma_dm[1][0], master_upgma_dm[1][2]])
# For combining on s4 and s5 (referencing original master index)
s3s4_s1s3 = np.mean([master_upgma_dm[0][3], master_upgma_dm[0][4], master_upgma_dm[2][3], master_upgma_dm[2][4]])
We can now fill in all of the distances in our iteration 2 distance matrix.
iter2_dm = [[0.0, s2_s1s3, s3s4_s1s3],
[s2_s1s3, 0.0, 5.5],
[s3s4_s1s3, 5.5, 0.0]]
iter2_dm = DistanceMatrix(iter2_dm, iter2_ids)
print(iter2_dm)
Step 4: There is still more than one value below the diagonal, so we start a new iteration by again repeating the process.
Iteration 3¶
Step 1: The smallest non-zero distance in the above matrix is now between (s1, s3) and s2. So, we'll draw that clade and set each branch length to half of the distance.
Step 2: We'll next create a new distance matrix where the clade (s1, s3) and the sequence s2 are now represented by a single clade, ((s1, s3), s2).
iter3_ids = ['((s1, s3), s2)', '(s4, s5)']
iter3_dm = [[None, None],
[None, 0.0]]
Step 3: We'll now fill in the values from the new clade to each of the existing sequences (or clades). This is now the mean of six distances. Why?
s1s2s3_s4s5 = np.mean([master_upgma_dm[0][3], master_upgma_dm[0][4],
master_upgma_dm[2][3], master_upgma_dm[2][4],
master_upgma_dm[1][3], master_upgma_dm[1][4]])
We fill this value into our iteration 3 distance matrix.
iter3_dm = [[0.0, s1s2s3_s4s5],
[s1s2s3_s4s5, 0.0]]
iter3_dm = DistanceMatrix(iter3_dm, iter3_ids)
print(iter3_dm)
Step 4: At this stage, there is only one distance below the diagonal in our distance matrix. So, we can use that distance to draw the final branch. This will connect our two deepest clades, ((s1, s3), s2) and (s4, s5), which will give us our final UPGMA tree.
SciPy contains an implementation of UPGMA that we can apply to our existing distance matrices, and we can then visualize the resulting trees with ete3. IAB provides a wrapper function that will give this an interface that is convenient to work with. If you'd like to see what the wrapper function is doing, using the psource IPython magic function as we have in other places in the text.
from iab.algorithms import tree_from_distance_matrix, guide_tree_from_sequences
%psource tree_from_distance_matrix
#%psource guide_tree_from_sequences
from scipy.cluster.hierarchy import average
help(average)
Let's compute and visualize UPGMA trees for the two distance matrices that we created above. How do these trees compare to one another? Does one look more or less correct than the other (they may or may not, depending on the random sample of sequences that are being compared).
One thing to be aware of as you start visualizing trees is that the vertical order (in the default TreeStyle being used here) doesn't have biological meaning, it's purely a visualization component. You can rotate the branches descending from any node in the tree freely.
kmer_tree = tree_from_distance_matrix(kmer_dm,
metric='upgma')
print(type(kmer_tree))
ete3.Tree(str(kmer_tree), format=1).render("%%inline", tree_style=ts)
jc_corrected_hamming_tree = tree_from_distance_matrix(jc_correct_dm(hamming_dm),
metric='upgma')
ete3.Tree(str(jc_corrected_hamming_tree), format=1).render("%%inline", tree_style=ts)
Which is better?
Which produces a higher percentage of agreement between leaf pairs?
How could we compute that?¶
As mentioned above, UPGMA has a rather complex sounding name: Unweighted Pair Group metric with Arithmetic mean.
- The Unweighted term indicates that all tip-to-tip distances contribute equally to each average that is computed (no weighted averages are being computed).
- The Pair Group term implies that all internal nodes, including the root node, will be strictly bifurcating, or descent to exactly two other nodes (either internal or terminal).
- The Arithmetic mean term implies that distances to each clade are the mean of distances to all members of that clade.
One invalid assumption that is made by UPGMA is inherent in Step 1, where each branch connecting the internal node to a tip is set to half of the length between the tips.¶
- This assumes the mutation rates are constant throughout the tree, or in other words that the tree is ultrametric.
- This is not likely to be the case in the real world, as different lineages in the tree might be undergoing different selective pressures, leading to different rates of evolution.
Neighboring joining is a distance-based phylogenetic reconstruction approach that does not assume ultrametricity.¶
- Neighbor-joining methods apply general data clustering techniques to sequence analysis using genetic distance as a clustering metric.
The simple neighbor-joining method produces unrooted trees, but it does not assume a constant rate of evolution (i.e., a molecular clock) across lineages.
http://scikit-bio.org/docs/latest/generated/skbio.tree.nj.html
- https://en.wikipedia.org/wiki/Neighbor_joining
- http://evolution.genetics.washington.edu/phylip/doc/neighbor.html
nj_tree = tree_from_distance_matrix(jc_correct_dm(hamming_dm),
metric='nj')
ete3.Tree(str(nj_tree), format=1).render("%%inline", tree_style=ts)
The following is a rooted tree, which means that it includes an assumption about the last common ancestor of all sequences represented in the tree.
https://en.wikipedia.org/wiki/Phylogenetic_tree#Rooted_tree
https://en.wikipedia.org/wiki/Tree_(graph_theory)#Rooted_tree
A rooted phylogenetic tree (see two graphics at top) is a directed tree with a unique node, the root, corresponding to the (usually imputed) most recent common ancestor of all the entities at the leaves of the tree.¶
- The root node does not have a parent node, but serves as the parent of all other nodes in the tree.
- The root is therefore a node of degree 2 while other internal nodes have a minimum degree of 3 (where "degree" here refers to the total number of incoming and outgoing edges).
The most common method for rooting trees is the use of an uncontroversial outgroup close enough to allow inference from trait data or molecular sequencing, but far enough to be a clear outgroup.
An unrooted tree, like the following, doesn't include an assumption about the last common ancestor of all sequences:¶
- Unrooted trees illustrate the relatedness of the leaf nodes without making assumptions about ancestry. T
- They do not require the ancestral root to be known or inferred.
- Unrooted trees can always be generated from rooted ones by simply omitting the root.
By contrast, inferring the root of an unrooted tree requires some means of identifying ancestry. This is normally done by including an outgroup in the input data so that the root is necessarily between the outgroup and the rest of the taxa in the tree, or by introducing additional assumptions about the relative rates of evolution on each branch, such as an application of the molecular clock hypothesis.
- Horizontal gene transfer (HGT) or lateral gene transfer (LGT)is the movement of genetic material between unicellular and/or multicellular organisms other than by the ("vertical") transmission of DNA from parent to offspring.
- Horizontal gene transfer is the primary mechanism for the spread of antibiotic resistance in bacteria, plays an important role in the evolution of bacteria that can degrade novel compounds such as human-created pesticides and in the evolution, maintenance, and transmission of virulence.
- Most thinking in genetics has focused upon vertical transfer, but horizontal gene transfer is important, and among single-celled organisms is perhaps the dominant form of genetic transfer.
Tree of life showing vertical and horizontal gene transfers:
A phylogenetic network or reticulation is any graph used to visualize evolutionary relationships (either abstractly or explicitly) between nucleotide sequences, genes, chromosomes, genomes, or species.¶
- They are employed when reticulation events such as hybridization, horizontal gene transfer, recombination, or gene duplication and loss are believed to be involved.
- They differ from phylogenetic trees by the explicit modeling of richly linked networks, by means of the addition of hybrid nodes (nodes with two parents) instead of only tree nodes (a hierarchy of nodes, each with only one parent).
- Phylogenetic trees are a subset of phylogenetic networks.
Phylogenetic network of HVS-I variation in mitochondrial DNA haplogroup R0 and R0b
Mini lab¶
Improving simulation assumptions
Q1¶
- Does our ability to re-construct improve with sequence length? Show a couple tests to see, using hamming/kmer and/or upgma/nj. Note: you don't need a totally computational metric of quality, but at least couple plots where you manually assess quality for different starting parent sequence lengths.
- Does it matter the method?
- If it changes, is the function linear, increasing, decreasing?
# Do coding here,
# snap screenshots of code/plot,
# and write a paragraph in Canvas,
# with screenshots embedded in the text entry box
# little helper function
def mean_distance(D):
sum = 0
for i in range(len(D)):
for j in range(len(D)):
if i != j:
sum += D[i][j]
return sum/(len(D)*(len(D)-1))
sequences = []
sequences_aligned = []
kmer_av_distances = []
hamming_av_distances = []
import pandas as pd
print(hamming_av_distances)
print(kmer_av_distances)
# entries go up by 50 sequences
def mean_line_graph(hamming_dists, kmer_dists):
cols = []
for i in range(len(hamming_dists)):
val = (i+1)*50
cols.append(str('%d seqs'%(val)))
df = pd.DataFrame({
'hamming': hamming_dists,
'kmer': kmer_dists
},index=cols)
df.plot.line()
def create_evolutions(seq_count, indel_prob):
# create a data set
sequence = random_sequence(skbio.DNA, seq_count)
sequences = evolve_generations(sequence, generations=3,
substitution_probability=0.1,
indel_probability=indel_prob,
increased_rate_probability=0.1,
verbose=False)
return sequences
def create_kmer_dm(sequences, distances, length=50):
kmer_dm = DistanceMatrix.from_iterable(sequences,
metric=kmer_distance,
key='id')
_title = '3mer distances between sequences length %d' %(length)
_ = kmer_dm.plot(cmap='Greens',
title=_title)
df = kmer_dm.to_data_frame()
df = df.values.tolist()
distances.append(mean_distance(df))
return kmer_dm
def create_hamming_dm(sequences, distances, length=50):
if indel_probability == 0:
sequences_aligned = sequences
else:
gpa = partial(global_pairwise_align_nucleotide,
penalize_terminal_gaps=True)
sequences_aligned = progressive_msa(sequences,
pairwise_aligner=gpa)
hamming_dm = DistanceMatrix.from_iterable(sequences_aligned,
metric=hamming,
key='id')
_title='Hamming distances between sequences length %d'%(length)
_ = hamming_dm.plot(cmap='Greens',
title=_title)
df = hamming_dm.to_data_frame()
df = df.values.tolist()
hamming_av_distances.append(mean_distance(df))
return hamming_dm
def create_kmer_tree(kmer_dm):
kmer_tree = tree_from_distance_matrix(jc_correct_dm(kmer_dm),
metric='nj')
ete3.Tree(str(kmer_tree), format=1).render("%%inline", tree_style=ts)
def create_hamming_tree(kmer_dm):
kmer_tree = tree_from_distance_matrix(jc_correct_dm(kmer_dm),
metric='nj')
ete3.Tree(str(kmer_tree), format=1).render("%%inline", tree_style=ts)
def create_both_trees(kmer,hamming):
create_kmer_tree(kmer)
create_hamming_tree(hamming)
def run_all_graphs(iter_counts, indel_prob, kmer_av, hamming_av):
for seq_count in iter_counts:
seqs = create_evolutions(seq_count, indel_prob)
hdm = create_hamming_dm(seqs, hamming_av, seq_count)
kdm = create_kmer_dm(seqs, kmer_av, seq_count)
create_both_trees(kdm,hdm)
hamming_av_distances = []
kmer_av_distances = []
seqs = []
for i in range(8):
seqs.append((i+1)*50)
run_all_graphs(seqs, 0, kmer_av_distances, hamming_av_distances)
mean_line_graph(hamming_av_distances, kmer_av_distances)
Q2¶
- Does our ability to re-construct improve with mutation rates that are uneven across the DNA segment, e.g., some regions being hyper-variable, and others being less so? No plots are needed for this one, just arm-chair thinking. Write a pagagraph in Canvas.